rm(list = ls())
set.seed(1111)
library(MatchIt)
library(dplyr)
library(ggplot2)
library(haven)
library(stargazer)
library(tableone)
library(Hmisc)
library(cem)
library(Matching)
library(car)
library(Zelig)
library(RColorBrewer)
library(car)
library(arm)
setwd("~/Dropbox/AttorneyMainData")
data <- read_dta("Attorney.dta")
data$hysLaw <- car::recode(data$SchoolLawyerArgue,"'Harvard'=1;'Yale'=1;'Stanford'=1;'Columbia'=1;'Chicago'=1;else=0")
data$demPres <- car::recode(data$SGOrals,"'McCree'=1;'Lee'=0;'(Acting) Fried'=0;'Fried'=0;'(Acting) Bryson'=0;'Starr'=0;'Days'=1;'(Acting) Dellinger'=1;'Waxman'=1;'(Acting) Underwood'=0; 'Olson'=0;'(Acting) Clement'=0;'Clement'=0;'Garre'=0;'(Acting) Kneedler'=1;'Kagan'=1;'(Acting) Katyal'=1;'Verrilli'=1;'(Acting) Gershengorn'=1;'(Acting) Francisco'=1;'(Acting) Wall'=1; 'Francisco'=1")
##Bryson was acting for both Bush 1 and Reagan
data$demPres[data$SGOrals=="(Acting) Bryson" & data$term=="1992"] <- 1
data$SGOralsCollapseF <- as.factor(data$SGOralsCollapse)
data$justiceNameF <- as.factor(data$justiceName)
data$ActingSGOrals <- as.numeric(factor(data$ActingSGOrals, ordered=F))-1
data$TopLawSchool <- as.numeric(factor(data$TopLawSchool, ordered=F))-1
data$DCFirm <- as.numeric(factor(data$DCFirm, ordered=F))-1
data$ClerkDummy <- as.numeric(factor(data$ClerkDummy, ordered=F))-1
data$OppPetOrResp <- as.numeric(factor(data$OppPetOrResp, ordered=F))-1
data$CriminalDummy <- as.numeric(factor(data$CriminalDummy, ordered=F))-1
data$CivLibDummy <- as.numeric(factor(data$CivLibDummy, ordered=F))-1
data$EconDummy <- as.numeric(factor(data$EconDummy, ordered=F))-1
data$Acela <- as.numeric(factor(data$Acela, ordered=F))-1
data$LawyerArgueOSG <- as.numeric(factor(data$LawyerArgueOSG, ordered=F))-1
data$GenderLawyerArgue <- as.numeric(factor(data$GenderLawyerArgue, ordered=F))-1
data$RaceLawyerArgue <- as.numeric(factor(data$RaceLawyerArgue, ordered=F))-1
data$IdeologyOppCounsel <- as.numeric(factor(data$IdeologyOppCounsel, ordered=F))-1
ctLevel <- subset(data, ct_level==1)
summary(data)
dim(data)

xvars<-c("SGOralsCollapseF","ActingSGOrals","hysLaw","DCFirm","ClerkDummy","IdeologyOppCounsel","Rehnquist","Roberts","OppPetOrResp","demPres","justiceNameF","MQCompatability","lcConservative","MQTerm")
data_nomiss <- data %>%  # MatchIt does not allow missing values
  dplyr::select(PriorOral, JVote, one_of(xvars)) %>%
  na.omit()
data_nomiss <- data.frame(data_nomiss)
xvars_case<-c("SGOralsCollapseF","ActingSGOrals","hysLaw","DCFirm","ClerkDummy","IdeologyOppCounsel","Rehnquist","Roberts","OppPetOrResp","demPres","MQCompatability","lcConservative")
data_case_nomiss <- ctLevel %>%  # MatchIt does not allow missing values
  dplyr::select(PriorOral, WhoWon, minWinning, unan, nearUnan, one_of(xvars_case)) %>%
  na.omit()
data_case_nomiss <- data.frame(data_case_nomiss)

match_justice_ext <- matchit(PriorOral ~ SGOralsCollapseF+ActingSGOrals+hysLaw+DCFirm+ClerkDummy+IdeologyOppCounsel+demPres+Rehnquist+Roberts+OppPetOrResp+justiceNameF,
                             method = "cem", data = data_nomiss)
summary(match_justice_ext)
dta_justice_ext <- match.data(match_justice_ext)

match_case_ext <- matchit(PriorOral ~ SGOralsCollapseF+ActingSGOrals+hysLaw+DCFirm+ClerkDummy+IdeologyOppCounsel+demPres+Rehnquist+Roberts+OppPetOrResp,
                          method = "cem", data = data_case_nomiss)
summary(match_case_ext)
dta_case_ext <- match.data(match_case_ext)

i1_justice_m_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres, family="binomial", data = data)
i1_justice_m_match_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
i1_justice_m <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp*IdeologyOppCounsel+demPres+justiceNameF, family="binomial", data = data)
i1_justice_m_match <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp*IdeologyOppCounsel+demPres+justiceNameF, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
i1_justice_m_noFE1 <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp*IdeologyOppCounsel+demPres, family="binomial", data = data)
i1_justice_m_match_noFE1 <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp*IdeologyOppCounsel+demPres, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
stargazer(i1_justice_m_noFE,i1_justice_m,i1_justice_m_match_noFE,i1_justice_m_match,i1_justice_m_noFE1,i1_justice_m_match_noFE1,type="text",digits=2,star.cutoffs = c(.05))

i2_justice_m_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+MQCompatability, family="binomial", data = data)
i2_justice_m_match_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+MQCompatability, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
i2_justice_m <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+MQCompatability+justiceNameF, family="binomial", data = data)
i2_justice_m_match <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+MQCompatability+justiceNameF, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
###SI, Table A6
stargazer(i2_justice_m_noFE,i2_justice_m,i2_justice_m_match_noFE,i2_justice_m_match,digits=2,star.cutoffs = c(.05))

i3_justice_m_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+lcConservative*MQTerm, family="binomial", data = data)
i3_justice_m_match_noFE <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+lcConservative*MQTerm, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
i3_justice_m <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+lcConservative*MQTerm+justiceNameF, family="binomial", data = data)
i3_justice_m_match <- glm(JVote~PriorOral+hysLaw+DCFirm+ClerkDummy+Rehnquist+Roberts+OppPetOrResp+IdeologyOppCounsel+demPres+lcConservative*MQTerm+justiceNameF, family="binomial", data = match.data(match_justice_ext), weights=match.data(match_justice_ext)$weights)
###SI, Table A7
stargazer(i3_justice_m_noFE,i3_justice_m,i3_justice_m_match_noFE,i3_justice_m_match,digits=2,star.cutoffs = c(.05))


##simulate
set.seed(1111)
dd <- data
aa <- sim(i2_justice_m_noFE)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.rm=T) + acoef[,"demPres"]*0  )
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.rm=T) + acoef[,"demPres"]*0 )
ideo1 <- treated-untreated
mean(ideo1)

dd <- data
aa <- sim(i2_justice_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.rm=T) + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.rm=T) + acoef[,"justiceNameFJPStevens"]*1)
ideo2 <- treated-untreated

dd <- match.data(match_justice_ext)
aa <- sim(i2_justice_m_match_noFE)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.omit=T))
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.omit=T))
ideo3 <- treated-untreated

dd <- match.data(match_justice_ext)
aa <- sim(i2_justice_m_match)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.omit=T) + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQCompatability"]*median(dd$MQCompatability,na.omit=T) + acoef[,"justiceNameFJPStevens"]*1)
ideo4 <- treated-untreated
mean(ideo4)

##mode of lcConservative is 0, which is why that and the interaction are not in the following lines of code.
dd <- data
aa <- sim(i3_justice_m_noFE)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1 + acoef[,"MQTerm"]*median(dd$MQTerm,na.rm=T) + acoef[,"demPres"]*0  )
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1 + acoef[,"MQTerm"]*median(dd$MQTerm,na.rm=T) + acoef[,"demPres"]*0 )
ideo5 <- treated-untreated
mean(ideo5)

dd <- data
aa <- sim(i3_justice_m)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQTerm"]*median(dd$MQTerm,na.rm=T) + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQTerm"]*median(dd$MQTerm,na.rm=T) + acoef[,"justiceNameFJPStevens"]*1)
ideo6 <- treated-untreated

dd <- match.data(match_justice_ext)
aa <- sim(i3_justice_m_match_noFE)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0  + acoef[,"MQTerm"]*median(dd$MQTerm,na.omit=T))
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQTerm"]*median(dd$MQTerm,na.omit=T))
ideo7 <- treated-untreated

dd <- match.data(match_justice_ext)
aa <- sim(i3_justice_m_match)
acoef <- coef(aa)
head(acoef)
treated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*1 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQTerm"]*median(dd$MQTerm,na.omit=T) + acoef[,"justiceNameFJPStevens"]*1)
untreated <- invlogit(acoef[,"(Intercept)"] + acoef[,"PriorOral"]*0 + acoef[,"hysLaw"]*median(dd$hysLaw,na.omit=T) + acoef[,"DCFirm"]*median(dd$DCFirm,na.omit=T) + acoef[,"ClerkDummy"]*median(dd$ClerkDummy,na.omit=T) + acoef[,"Rehnquist"]*1 + acoef[,"Roberts"]*0 + acoef[,"OppPetOrResp"]*median(dd$OppPetOrResp,na.omit=T) + acoef[,"IdeologyOppCounsel"]*1  + acoef[,"demPres"]*0 + acoef[,"MQTerm"]*median(dd$MQTerm,na.omit=T) + acoef[,"justiceNameFJPStevens"]*1)
ideo8 <- treated-untreated
mean(ideo8)

###Manuscript, Figure 6
pdf("Ideology.pdf",width=15,height=6.6,paper='special') 
par(mar=c(5.1,2, 4.1, 2.1),mfrow=c(1,2),oma=c(0,7,0,0))
plot(NA,xlim=c(-.05,.2),ylim=c(.5,2.5),ylab="",xlab="ATT", main="Ideological Congruity",yaxt="n")
axis(2, at=c(2,1),labels=c("No Fixed Effects","Justice Fixed Effects"), las=2)
points(y=c(2.1,1.1),x=c(mean(ideo1),mean(ideo2)), pch=19,cex=1.5,col="gray70")
points(y=c(1.9,.9),x=c(mean(ideo3),mean(ideo4)), pch=18,cex=1.75)
segments(y0=2.1,y1=2.1,x0=quantile(ideo1,.025),x1=quantile(ideo1,.975),lwd=2,col="gray70")
segments(y0=1.1,y1=1.1,x0=quantile(ideo2,.025),x1=quantile(ideo2,.975),lwd=2,col="gray70")
segments(y0=1.9,y1=1.9,x0=quantile(ideo3,.025),x1=quantile(ideo3,.975),lwd=2)
segments(y0=.9,y1=.9,x0=quantile(ideo4,.025),x1=quantile(ideo4,.975),lwd=2)
abline(v=0,lwd=3,col="gray60")

plot(NA,xlim=c(-.05,.2),ylim=c(.5,2.5),ylab="",xlab="ATT", main="Interaction Term",yaxt="n")
points(y=c(2.1,1.1),x=c(mean(ideo5),mean(ideo6)), pch=19,cex=1.5,col="gray70")
points(y=c(1.9,.9),x=c(mean(ideo7),mean(ideo8)), pch=18,cex=1.75)
segments(y0=2.1,y1=2.1,x0=quantile(ideo5,.025),x1=quantile(ideo5,.975),lwd=2,col="gray70")
segments(y0=1.1,y1=1.1,x0=quantile(ideo6,.025),x1=quantile(ideo6,.975),lwd=2,col="gray70")
segments(y0=1.9,y1=1.9,x0=quantile(ideo7,.025),x1=quantile(ideo7,.975),lwd=2)
segments(y0=.9,y1=.9,x0=quantile(ideo8,.025),x1=quantile(ideo8,.975),lwd=2)
abline(v=0,lwd=3,col="gray60")
legend("bottomright",c("Full Dataset","Matched Dataset"),       col=c("gray70","black"),pch=c(19,18),pt.cex=c(1,1.3))
par(mar=c(5.1,4.1, 4.1, 2.1),oma=c(0,0,0,0))
par(mfrow=c(1,1))
dev.off()
